Time evolution of simple molecules during proto-star collapse 
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Abstract 

We study the formation and evolution of several molecules in a collapsing interstel- 
lar cloud using a reasonably large reaction network containing more then four hundred 
atomic and molecular species. We employ a time dependent, spherically symmetric, hy- 
drodynamics code to follow the hydrodynamic and chemical evolution of the collapsing 
cloud. The flow is assumed to be self-gravitating. We use two models to study the hy- 
drodynamic evolution: in the first model, we inject matter into an initially low density 
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region and in the second model, we start with a constant density cloud and let it collapse 
due to self-gravity. We study the evolution of the central core for both the cases. We 
include the grain chemistry to compute the formation of molecular hydrogen and carried 
out the effect of gas and grain chemistry at each time step. We follow the collapse for more 
than 10 14 s (about 3 million years) and present the time evolution of the globally averaged 
abundances of various simple but biologically important molecules, such as glycine, ala- 
nine etc. We compare our results with those obtained from observations found that for 
lighter molecules the agreement is generally very good. For complex molecules we tend 
to under predict the abundances. This indicates that other pathways could be present to 
form these molecules or more accurate reaction rates were needed. 
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1. Introduction 

More than 125 species of molecules have been observed in the interstellar clouds and 
star forming regions. Among them, over half are organic. Serious efforts have been made 
over the years to investigate the formation of such molecules in cool interstellar clouds 
in frigid conditions (Hasegawa et al., 1992; Hasegawa and Herbst, 1993; Leung et al., 
1984; Prasad and Huntress, 1980a, 1980b). It is now quite certain that the most impor- 
tant building block, namely, the molecular hydrogen (H 2 ) and some of the other lighter 
molecules must be produced in the presence of grains (Gould and Salpeter, 1963; Hollen- 
bach and Salpeter, 1971; Hollenbach et al., 1971). Several analytical and numerical works 
have successfully shown how the molecular hydrogen may have been produced (Biham 
et al., 2001). A number of results are present in the literature where hydrodynamic and 
chemical evolutions have been attempted simultaneously. For example, Shalabiea and 
Greenberg (1995) used the pseudo as well as partially real time-dependent models for 
the hydrodynamical evolution. In the pseudo time-dependent method, they assumed a 
constant density and temperature of the cloud using which the chemical evolution was 
computed. In their time-dependent model, they included the density and temperature 
variations throughout the cloud. However, in their initial approach to the time-dependent 
modelling, they assumed a constant temperature but allowed only the density to vary. 
Ceccarelli et al. (1996) used the "inside-out", isothermal, spherical collapse model of Shu 
(1977) and coupled it with a time-dependent chemical evolution code. They included the 
heating and the cooling processes with an emphasis on the line emission. Shematovich et 
al. (1997) used Zeus 2D code which included the heating and the cooling. They present 
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ID hydrodynamic and chemo-dynamical evolution of the proto-stellar cloud illuminated 
by the diffused interstellar UV radiation. They solved the equations of chemical kinetics, 
hydrodynamics and thermal balance simultaneously. In Lim et al. (1999) 2D numerical 
code was developed using the adaptive grid technique. Here, 454 reactions among 42 
atomic and molecular chemical species were taken including the basic elements like H, 
He, C, N, O and a representative low ionization potential metal Na. At each grid point, 
the chemical evolution was followed by a calculation of the reaction rates using the lo- 
cal conditions obtained from the hydrodynamical flow. They primarily concentrated on 
the diffused clouds and emphasized the interfaces of the interstellar media and resulting 
dynamical mixing. Aikawa et al. (2005) studied time-dependent evolution of Bonner- 
Ebert spheres by assuming clouds having a specific parameter a which is the ratio of 
the gravitational force to the pressure force. Recently, Acharyya et al. (2005) solved 
the Master equations and rate equations of Biham et al. (2001) for various cloud pa- 
rameters and followed the evolution of H 2 as a function of time. Both this work and 
the earlier works of Chakrabarti and Chakrabarti (2000a) employed steady state matter 
distribution and assumed that the density and the temperature distributions at a given 
radial distance do not change with time. Chakrabarti and Chakrabarti (2000a) used a 
large number of species and the reaction rates were taken from the UMIST data base. 
Some of the reaction rates which were not available in the literature were assumed to 
be similar to other two body reactions. Subsequently, these new and assumed reaction 
rates were parametrized (with reaction rates up to a thousand times smaller compared 
to Chakrabarti and Chakrabarti (2000a) to include the effect of the size of the reactant 
molecules (Chakrabarti and Chakrabarti, 2000b). It was shown that even under frigid 
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and tenuous conditions of the interstellar media, a significant and perhaps a detectable 
amount of simple amino acids and even important ingredients of DNA molecule (such 
as adenine) may form. Ceccarelli et al. (2000) estimated the upper limit of the abun- 
dance of glycine to be about 10~ 10 (cooler outer cloud) to 7 x 10~ 9 (hot core). Kuan et 
al. (2003) estimated the fractional abundance of glycine to be 2.1 x 10~ 10 for Sgr B2, 
1.5 x 10~ 9 for Orion, and 2.1 x 10~ 10 for W51. These numbers are comparable to what 
was predicted in Chakrabarti and Chakrabarti (2000a), however, there are clearly some 
debate on the possible pathways for the formation of glycine with the route followed in 
the Chakrabarti and Chakrabarti (2000a, 2000b). Similarly, there are also some debate 
on whether glycine is actually observed (Hollis et al., 2003; Snyder et al., 2005). A few 
other relevant results which may be mentioned in passing are as follows. Tarafdar et 
al. (1985) presented a model of chemical and dynamical evolution of isolated, initially 
diffused and quiescent interstellar clouds. A semi-empirically derived dependence of the 
observed cloud temperatures on the visual extinction and density was used in this work. 
Sorrell (2001) outlined a theoretical model for the formation of the interstellar amino acids 
and sugars. In this model, first ultraviolet photolysis creates a high concentration of free 
radicals in the mantles and the heat input due to the grain-grain collision causes radicals 
to react chemically with another to build complex organic molecules. Bernstein et al. 
(2002) reported a laboratory demonstration that the complex bio-molecules like glycine, 
alanine etc. are naturally formed from the ultraviolet photolysis of the interstellar grains. 
Munoz Caro et al. (2002) report the detection of amino acids even at room temperature 
on an interstellar ice analogue that was irradiated with ultraviolet light in a high vacuum 
at 12K. Altogether sixteen amino acids were identified. The results demonstrate that a 
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spontaneous generation of amino acids in the interstellar medium is possible, supporting 
the suggestion that prebiotic molecules could have been delivered to the early earth by 
cometary dust, meteorites or interplanetary dust particles. 

In this backdrop of this observational and experimental status, we carry out our inves- 
tigation by improving earlier work by first incorporating the accurate grain chemistry as 
elaborated in Acharyya et al. (2005) and then by actually combining the results of a time 
dependent hydrodynamics code with the chemical evolution code to see how the abun- 
dances vary with the grid locations. We also chose initial conditions very much different 
from the earlier studies. We used two different but realistic models. In one model (Model 
A), the cloud matter is injected through the grid boundary at a constant speed and the 
cloud as well as the core are allowed to form ab initio. In the other model (Model B), the 
computational grid area is chosen to be the central part of a much larger spherical cloud 
of constant density and temperature. The rate at which the interstellar matter enters into 
the grid (i.e., the rate at which the larger cloud is evacuated) depends on the gravitational 
pull between the inner cloud within computatioal grid and the outer cloud outside of the 
computational grid. In following the chemical evolution, we used the available standard 
reaction rates for most of the reactions, but the rates of some of the very complex molecule 
formation are still very much uncertain and as such we present the evolution of the mass 
fractions only for simpler bio-molecules. In the next section, we present the hydrody- 
namic equations which govern the cloud collapse. In §3, we present the hydrodynamical 
and chemical evolution. In §4, we discuss in detail the nature of the cloud models which 
we simulate and how the chemical evolution code is used in conjunction with the results 
of the hydrodynamic simulation. In §5, we present the results. Finally, in §6, we make 
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2. Hydrodynamic Equations Governing the Cloud Collapse 

The time dependent simulations have been carried out by various authors using both 
the Lagrangian and Eulerian methods. For instance, Prasad et al. (1991) and Tarafdar 
et al. (1985) used a Lagrangian scheme with a semi-empirical approach for the energy 
equation. Shalabiea and Greenberg (1995) used a partially time dependent model in 
which only the densities are allowed to be time dependent during collapse of the cloud. 
We have used a finite difference Eulerian scheme (upwind scheme) in which we difference 
the one dimensional Euler equations along the radial grid. We consider a self-gravitating, 
collapsing, spherically symmetric gas cloud. We choose the co-ordinate system to be 
(r, 9, <p) with origin at the center of the proto-star. Since we would be interested in 
the spherical case, we ignore the 6 and <fi component motions and concentrate only on 
the radial equation in this paper. The Eulerian equations of hydrodynamics written in 
spherical coordinates are given by: 




(1) 



and 



dpv r 1 d(pv 2 r 2 ) 




<9$ dp 



(2) 



dt r 2 dr 



Here, $ is the gravitational potential [= —GM(r)/r] where, M(r) is the mass of the cloud 
inside radial distance r. 
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We assume the ideal gas equation to be, 

p = pkT//im p , (3) 

where, k is the Boltzmann constant, T is the local temperature, \i is the mean molecular 
weight, m p is the proton mass. 

3. Hydro dynamical and Chemical Evolution 

The equations presented in the previous Section have been solved to study full time- 
dependent behaviour of the spherical flow. The actual solution depends on the model 
flow. In the present paper we used two types of initial conditions for the evolution of two 
hydrodynamic models. These are coupled with chemical evolution code to obtain time 
dependent chemical composition. 

3.1. Hydrodynamical Models 

(a) Model A: Here we start with a grid of size r out . We assume that initially the cloud 
contains a negligible amount of mass. At the outer boundary, we inject matter at a 
constant rate of ^jf = 4:%p out v ou trl ut , where, p out is the injection density and v out is the 
injection velocity at the outer boundary (r = r out ). Using this model, we mimic the 
formation of the cloud itself from a supply of matter from a large reservoir of diffused gas. 

The inner boundary is chosen at r = r in . The distance r out — r in is divided into 
N = 100 logarithmically equal spaced grid points. At each time step, the mass of the 
core is dynamically updated by the amount of matter that is getting inside r in . The rate 
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dt 



where, p in is the density at the inner grid point, v rin is the velocity at the inner grid point, 
r in is the radial distance of the inner grid point from the center of the molecular cloud. 
Similarly, the mass of the cloud is also dynamically updated. The flow is assumed to be 
isothermal throughout and T = 10K is assumed for concreteness. 

(b) Model B: In this model, we assume that we have a finite sized r = r res reservoir 
of matter and we are considering only the inner region of size r out for computational 
purpose. The matter is coming from the reservoir and it is evacuated to form the star itself. 
The inward velocity is computed self-consistently from the gravitational pull between the 
matter inside r out and that outside of r out . The whole cloud is assumed to be isothermal 
(with T = 10K) and has a constant density p int throughout. The size of the whole cloud 
is r res » r out . We assume that M total = M core + M doud + M out , where, M core is the core 
mass within r = r in which will increase at the same rate as in eq. (4), M doud is the mass of 
the cloud within the computational grid, i.e., in between r in and r out and M out is the mass 
in between the radius r res and r out . Thus, M out is the reservoir mass which is depleted as 
matter is injected within r = r out . The injection velocity is dynamically calculated from 
the following way. First, we note that an approximate expression for the acceleration of 
matter at the outer edge of the computational grid could be taken as, 



out 



G (Mdoud + M core ) 

Tri? 



(5a) 
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where, r m = 



(j' Tes + r out ) / 2. Hence the injection velocity is, 




(56) 



and the average density of the injected matter is, 



where, V out = 4/3n(r^ es — r% ut ). By this process, the large cloud would be gradually 
evacuated as the star is formed. The cloud mass M doud inside r in <r < r out will increase 
at the rate, 



where, v in and p in are respectively the velocity and density at the inner boundary r in . 

3.2. Chemical Evolution 

The chemical evolution code is similar to that used in Chakrabarti and Chakrabarti 
(2000a, 2000b) except in one very important way, namely, the way we handle the grain 
chemistry. In Chakrabarti and Chakrabarti (2000a, 2000b) a constant rate of recombina- 
tion of H + H — > H 2 was chosen following the prescription of Miller et al. (1997). In our 
present work, however, we incorporate the grain chemistry properly for the production of 
H 2 at each radius. Here, we followed the techniques used in Acharyya et al. (2005). We 
used only the olivine grains for simplicity. We also incorporate the grain size distribution 
according to Weingartner and Draine (2001a, 2001b) and sub-divided the grains into three 
major types (a) 5 A, (b) 75 A, and (c) 0.2/x so that the calculation of H 2 becomes faster. 



dM c i oud 
dt 



4:7T(p out v out rl ut - p in v in rl 



(6) 
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Whether or not the Master equation or the rate equation would be used depends on the 
flux of the infalling hydrogen (Acharyya et al., 2005) on the grain surface. We take the 
reaction rates from Miller et al. (1997) which supplies the rate constant k for a two body 
reaction as: 

k = a(T/300)^exp(-7/T)cm 3 s _1 , (7) 

where, a, (3, 7 are constants. The reactions for which the rates were not available in the 
UMIST data base, we take a conservative value of a = 10~ 10 cm~ 3 s _1 , (3 = 7 = like any 
other typical two body (open shell) reactions even though we may be using neutral-neutral 
reactions. The chemistry of bio-molecule formation in space is very much unknown and it 
is difficult to quantify the rate with any certainty. Thus the results obtained using these 
rates may have some errors. 

Two models have been run with different initial abundances of carbon. Generally, we 
keep the initial composition as was kept in Miller et al. (1997) (after converting to mass 
fractions as in Chakrabarti and Chakrabarti, 2000a), i.e., H:He:C:N:0:Na:Mg:Si:P:S:Cl:Fe 
= 0.64:0.35897:5.6 x 10 4 :1.9 x 10~ 4 :1.81 x 10~ 3 :2.96 x 10~ 8 : 4.63 x 10^:5.4 x 10~ 8 :5.79 x 
10- 8 :4.12 x 9~ 7 :9 x 10- 8 :1.08 x lO" 8 . This Set is our 'Set 1' initial abundance. In 'Set 2' 
initial abundance, we used a higher amount of carbon to simulate the chemical evolution 
inside a collapsing cloud in the environment of evolved stars. 

While we include all the reactions given in the UMIST data base, we wish to remind 
that since the interstellar medium is cooler, the exothermic reactions should be more 
favourable. The following major types of the reactions are included (Tielens, 2005): 
(i) Photochemistry: 
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The UV photons present in the diffuse ISM are a dominant destruction agent for small 
molecules. Inside a cloud, the radiation field will be attenuated by the dust. In the 
unshielded radiation field this type of reactions have a typical rate of 1CT 9 s _1 . These 
reactions mainly occur in the outer edge of the molecular cloud where UV photons are 
profuse in number. Deep inside the cloud some UV photons can be created due to radiative 
association, which could also lead to this type of reactions. The rate is calculated by, 

n pd = aexp[-bA v ], 

where, A v is the visual extinction due to the dust, a is the unshielded rate and b is the 
self-shielded rate. An example of this kind of reaction is, 

CH + hv -> C + H. 

(ii) Neutral-neutral: This type of reactions often possesses appreciable activation barri- 
ers because of the necessary bond breaking associated with the molecular re- arrangement. 
This type of reactions is important when the gas is warm, e.g., in stellar ejecta, in hot cores 
associated with proto-stars, in dense photo-dissociation regions associated with luminous 
stars, or in the post-shock regions. The rates are generally around 4 x 1CT 11 cm 3 s _1 
when no activation barrier is present. Reaction rates are calculated as, 

k = a(T/300fexp(--f/KT), 



where, a, (3, 7 are used to calculate the rate coefficients. One example of this type of 
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reaction is, 

H 2 + ^ OH + H. 

In general, the only neutral-neutral reactions that occur in the cold conditions of the dark 
clouds are those involving atoms or radicals, often with non-singlet electronic ground 
states that do not have activation barriers. One example of this type of reaction is, 

C + OH — > CO + H. 

(iii) Ion-molecule reactions: Exothermic ion-molecule reactions occur rapidly because 
the strong polarization-induced interaction potential can be used to overcome any acti- 
vation barrier energy involved. In the exothermic direction, this kind of reactions have a 
typical rate of 2 x 10~ 9 cm 3 s" 1 and the reaction rates are in the form 

k — a. 

One example of this kind of reaction is, 

H% + Hi — > H% + H. 

(iv) Charge transfer reactions: This type of reactions is of great importance for 
setting the ionization balance in the HII region. The charge exchange between O and H + 
is a very important reaction in the ISM because it ionizes the oxygen which is then able 
to participate in the chemistry of the ISM. This type of reactions has a typical rate of 
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10 9 cm 3 s 1 in the exothermic direction. An example is: 

H + + O^H + + . 

(v) Radiative association reactions: In this type of reaction, the product after the 
collision of two species is stabilized through the emission of a photon. The reactions of 
this type have highly reaction specific rates. The reaction rates are given by, 

k — a. 

One example of this type of reaction is, 

H + C -> CH + hv. 

(vi) Dissociative electron recombination reactions: This type of reactions involves 
in the capture of an electron by an ion to form a neutral in an excited electronic state 
that can dissociate. Typical rate of such a reaction is around 1CT 7 cm 3 s _1 and the rate 
coefficient is calculated by, 

K = a(T/300) p . 
One example of this kind of reaction is, 

OH + + e^O + H. 

(vii) Associative detachment reactions: In this case, an anion and an atom collide 
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and the neutral product stabilizes through an electron emission. This type of reactions 
has a rate of around 10^ 9 cm 3 s^ 1 in the exothermic direction. 



H + H- -> H 2 + e. 



(viii) Collisional association: In laboratory settings, three-body reactions generally 
dominates chemistry, 

A + B + M^AB + M 

with rates ~ 10~ 32 cm 6 s" 1 . These reactions generally have very little importance in the 
astrophysical environment except for dense gas near stellar photosphere or in dense(~ 
10 n cm~ 3 ) circumstellar disks. 

We do not give the reaction pathways since there would be too many of them. Gener- 
ally, pathways given in Chakrabarti and Chakrabarti (2000a, 2000b) have been used for 
synthesis of those complex molecules which are absent in the UMIST data base. 

4. Results 

Armed with a working time-dependent hydrodynamical evolution code and a chemical 
evolution code, we ran two models (Models A and B) of the interstellar cloud collapse. In 
both the Models we choose, r in = 3.98 x 10 13 cm and use 100 grid points which are spaced 
equally in the logarithmic scale along the radial direction. We assume that anything 
going inside r in will increase the mass of core while the the radius of the core remains 
r in . In Model A, we chose r out = 3.54 x 10 18 cm and we ran for both the Sets 1-2 as the 
initial compositions. We started our simulation with p int = 10~ 27 gm cm -3 in all the grids 
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(except on the boundary, see below). Thus, we start with a core of a negligible mass of 
(4/3)irrf n p int = 2.64 x 10 14 gm. We inject matter at v int = 50cm s -1 at each grid point in 
the outer boundary. 

In Model Bl, we combine the hydrodynamical code and the chemical evolution code 
using only the Set 1 as the initial composition. Here, we choose r res = 10 20 cm, r out = 
3.16 x 10 18 cm and r in = 3.98 x 10 13 cm and use equal logarithmic radial spacings as 
before. At the start of the simulation p int = 1CT 22 gm cm~ 3 and v int = 10 cm s _1 at each 
grid point. Initially, the core mass was (A/3)7rrf n p int = 2.64 x 10 19 gm in this Model. In 
Model B2, we choose r res = 3.16 x 10 18 cm, r out = 3.16 x 10 1? cm and other conditions are 
same as above. 

As far as the chemical evolution is concerned, we use a large network of about 4000 
reactions and solve 422 equations simultaneously each of which is to update the mass 
fraction of one single specie. However, in order to simplify the chemical evolution and to 
save the computational time, we do not run the chemical evolution code at each grid at 
each time step. Instead, we divide the entire region under consideration in ten logarith- 
mically equal spaced zones along the radial direction. Appropriately weighted averaged 
density is used in each such zone. The initial abundance of matter in each zone is updated 
using the existing abundance in that zone plus the abundance of the upstream matter 
advecting in and minus the abundance of the downstream matter advecting out. The 
time step in chemical evolution is totally dictated by the fastest reaction in the network. 
The evolution process is continued till the end of the hydrodynamic simulation. While 
we have time evolution of the abundance of each species at each zone, we plot only the 
global average of the abundances just to give an idea of how the average abundance is 
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evolving. 

4.1. Results of Models A and B with Set 1 as initial composition 

In this Model, we choose p out = lCT 22 gm cm 3 , v out = 10 4 cm s _1 . Thus, the rate of 
injected matter is given by M = 4:iip out v ou tr^ ut = 1.57 x 10 20 gm s -1 . This matter moves 
in due to self-gravity of the cloud and the attraction of the core. The simulation was 
carried out till a few times 10 6 yr. Figure la shows how the velocity profile of the flow 
changes with time. Time (in years) is marked on each curve. A flow which began with a 
constant velocity, eventually assumes an almost steady state (v(r) ~ r~ a with a ~ 4/5). 
The evolution of the density distribution is shown in Figure lb. Here too, an initially 
constant density distribution assumes a power-law distribution of p(r) ~ r - 3 / 2 toward 
the end of the simulation. Time (in years) is marked on each curve. Note that once the 
core becomes massive, it starts to evacuate the grids and thus the density over the entire 
cloud gradually decreases as is shown by the dot-dashed curve drawn at ~ 6.33 x 10 6 yr. 
After about t = 10 13 s, i.e., about 3 x 10 5 yr, densities on each grid started increasing after 
the empty grids are filled in by the inflowing matter. This is shown in Figure lc where 
the time evolutions of the densities on the 50th grid (solid curve) and the innermost grid 
(dotted curve) are shown. Towards the end, we note that the densities at these grids 
gradually decrease as the cloud starts to become empty. Meanwhile, the core starts to 
grow because of the mass accretion during this period. Figure Id shows how the mass of 
the core is evolving. 

In Model Bl, the velocity assumes almost steady state after some initial transient time 
as before. The evolution of velocity is shown in Figure 2a. In Figure 2b, the evolution of 
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Figure 1. Evolutions of (a) the velocity and (b) the density distribution inside the col- 
lapsing cloud of a Model A simulation. Also shown are the evolutions of the density at 
(c) the middle (solid) grid and the innermost grid (dotted) and (d) the growth of the 
proto-stellar core with time. The densities in the grid have a sharp rise after the initially 
empty grid is filled in by the infalling matter at about 10 13 s. From this time onward, the 
core also starts growing rapidly. 
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Figure 2. Same as Figs.l(a-d) except that Model Bl is used for simulation. In this case, 
the depletion of matter from a cloud of fixed mass takes place and the core mass of the 
cloud goes up initially rapidly and slowly afterwards. 
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density is shown. Initially, the density of the cloud increases and new matter is injected 
from the cloud external to the grid, but as the time proceeds, the core grows while the 
external cloud is evacuated. This results in eventual decrease in density of the cloud as 
shown by the dot-dashed curve in Figure 2b. Similarly, from Figure 2c, it is clear that the 
density at the middle grid (solid curve) and the innermost grid (dotted curve) will also 
go down due to the growth of the core mass. In Figure 2d, the evolution of core is shown. 
Unlike in the previous model, the core mass starts growing from the very beginning of 
the simulation, but the rate of growth is slowed down when the cloud mass is depleted. 

Our Model A and Model B are two complimentary models. In Model A, the initial 
grid was almost empty and thus the core mass was growing slowly. Eventually the matter 
injected at the our boundary caught up and the core grew rapidly. In Model B, the initial 
grid was full. After it is evacuated, matter from the reservoir is accreted slowly. Thus, 
initially, the core grows rapidly but the rate of growth slowed down in the late phase. 

We now show the variation of different species for Model A hydrodynamic simulation. 
Figures 3(a-d) show the time variation of the average abundances of several species using 
the initial composition as in Miller et al. (1997) (Set 1). In Figure 3a we show the 
evolution of H and H 2 . The mass fraction of the atomic hydrogen goes down since it is 
utilized to form H 2 and other hydrogenated species. Since the atomic hydrogen is very 
much reactive agent it goes to the molecular form easily and thus towards the end of our 
simulation almost all the H was found primarily in the form of H 2 . In Figure 3b, we 
show the abundances of 2 , H 2 and CO whose final mean mass fractions are around 
10~ 5 , 10~ 7 and 10~ 3 respectively. The production rate suddenly increases after the grids 
are filled up with the collapsing shell of matter and the densities go up. In Figure 3c, 
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Figure 3. Time variation of the spatially averaged mass fractions of (a) H and H 2 , (b) 
2 , H 2 and CO (c) HCN and iV# 3 and (d) C 2 H 5 OH, CH 3 CHO and CH 3 OH. Here, 
Model A cloud is chosen and the initial mass fractions are the same as that of Set 1. 
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Figure 4. Time variation of the spatially averaged mass fractions of two simple bio- 
molecules, such as glycine and alanine. The hydrodynamic and chemical model are the 
same as in the previous Figure. 



the time dependences of the abundances of HCN and NH 3 are shown. The final mass 
fractions of these species are around ~ 1CT 7 . In Figure 3d, we show the variations of 
C 2 H 5 OH, CH 3 CHO and CH z OH whose final mass fractions are ~ KT 15 , ~ KT 10 and 
~ 10~ 9 respectively. Finally, in Figure 4, we show the evolutions of some of the simple 
bio-molecules, such as glycine and alanine whose final mass fractions are found to be 
around ~ 1CT 13 to ~ 10~ 15 respectively. 

4.2. Collapse of clouds with other initial carbon abundances and comparison 
with observations 

Since the initial composition of the clouds is not known very accurately, we ran Model 
A for one more initial chemical abundance (Set 2), where we use the mass fraction of 
carbon to be 0.001, roughly twice as much as we used in Set 1 above. This may mimic the 
environment of an evolved star. We adjusted the mass fraction of helium accordingly so 
as to keep the sum of the mass fractions to be unity. As expected, the carbon containing 
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Figure 5. Time variation of the spatially averaged mass fractions of H 2 and CO when 
the initial carbon abundance is twice that in Set 1. 
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Figure 6. Time variation of the spatially averaged mass fractions of (a)H 2 and CO 
(b)C 2 H 5 OH, CH 3 CHO and CH 3 OH. Here Model Bl was chosen. 



species become more abundant in simulations with Set 2 as the initial abundance. In 
Figure 5, we present the variations of 2 , H 2 and CO with time for these two models. 
These are to be compared with those presented in Figure 3b. 

Figures 6(a-b) show the time variation of the average abundances of several species 
using the initial composition of Model Bl. In this connection, we wish to note that 
though we have been using the so-called freeze-out of the chemical species, the effect is 
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Time (sec) Time (sec) 



Figure 7. Evolution of the mass fraction of CO (a) without and (b) with the freeze-out 
effect taken into account. The results in the 1st, 4th, 7th and 10th shells are displayed. 
Note that in the 1st and the 4th shell, the freeze-out effect is prominent. However, 
the global averages in both the cases do not differ by any large margin since they are 
dominated by the outer shells. Evolution in (b) only up to the freeze-out time is depicted 
to show detailed variation. 

not very easy to see in our abundance plots. This is because we are plotting the globally 
averaged abundances which is dominated by very slowly evolving matter in the outermost 
cell. By the time the matter is dense and more evolved, it is compressed to a smaller 
volume and does not contribute much in the weighted average. However, when we plot 
the evolution for each cell, we see that the freeze-out effect can be seen in the innermost 
four cells. In Figures 7(a-b) we show the evolutions of the mass fraction of CO in the 
1st (innermost), 4th, 7th and 10th (outermost) shells of matter in the cloud when (a) we 
do not consider the freeze-out effect and when (b) we do consider the freeze-out effect. 
When no freeze-out is considered (Figure 7a), all the shells eventually reach a saturation 
value of the mass-fraction. This is reached when the enhancement of the CO through 
combination of C and O matches its reduction through hydrogenation. Different shells 
take different times depending on the average density of each shell. When we assume a 
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Table 1 

Comparison of our Results with the observed abundances 



Molecules 


Observed 


Model A 




Model B 


References 


a HCN 


6 X 10- 9 


1.4 x 10~ 9 


2.4 x 10 


-!) 


1.9 x 10" a 


Irvine and Hjalmarson (1983) 


a NH z 


5 x 10~ 8 


4.5 x 10" 9 


1.2 x 10 


-9 


3.9 x 10" 9 


Tollc ct al. (1981) 


a CO 


3 x itr 5 


1.8 x 10~ 5 


6.6 x 10 


-5 


2.6 x 10" s 


Allen and Knapp (1978) 


a H 2 


7 x icr 8 


4 x 10" 8 


2.2 x 10 


-8 


1.6 x 10- 8 


Snell ct al. (2000) 


a o 2 


itr 6 


3.3 x lO" 7 


1.9 x 10 


-7 


2.7 x 10- 7 


Gold et al. (2000) 


a CH 3 OH 


5 x 10- 10 


2.5 x lO" 11 


4.4 x 10~ 


11 


1.6 x 10- 11 


Priberg et al. (1988) 


a CH 3 CHO 


3 x icr 10 


1.3 x 10" 13 


9.1 x 10- 


13 


1.9 x 10- 13 


Matthews et al. (1985) 


b C2H 5 OH 


1.5 x 10~ 9 


5.8 x 10~ 17 


3.6 x 10- 


Ifi 


1.1 x 10~ 15 


Ohishi ct al. (1995) 


Alanine 




2.3 x lO" 17 


8.3 x 10- 


17 


6 x 10- 17 




Glycine 


10 -10 c 


1.7 x 10- 14 


2.9 x 10- 


14 


1.7 x 10- 14 


Ceccarelli et al. (2000) 




7 x lO" 9 d 










and 




(0.21 - 1.5) x 10~ 9 










Kuan ct al. (2003) 



a Observed near TMC-1 (10K) 

b Observed near Orion KL (70K) 

c upper limit (cold cloud) ; d upper limit (hot core) 



freeze-out, the mass fractions in the innermost four shells start to decrease after reaching 
the saturation while the other shells do not decrease within the evolution time of the 
cloud. The evolution of the global average follows initially the evolution of the 10th shell, 
but afterwards, from around t = 2 x 10 13 s, the deviation occurs as the global average is 
dominated by denser intermediate shells. 

In Table 1, we present the results of Model A with Sets 1 and 2 chemical compositions, 
and the result from Model B using Set 1 composition (i.e., Model Bl). We also compare 
the observed abundances of some of the species. For comparison, we converted mass 
fraction of a species into ratio of the number density of that species and that of H 2 so 
that our results are directly comparable with observables. The references from where the 
observed abundances have been obtained are provided in Column 6. 

From the Table, it is clear that our results in Sets 1-2 are generally close to the ob- 
servational results for most of the lighter species. Since all the C and O are primarily 
channeled into the formation of CO, production of other and more complex hydrocar- 
bons by our method is not very efficient. For instance, for C 2 H 5 OH, glycine, alanine, 
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abundances obtained by our simulations are very low compared to the observed or 'upper 
limits' obtained from observations. According to Ohishi et al. (1995), C2H5OH is difficult 
to form in the gas phase. This is thought to be synthesized on the grain surfaces and 
then evaporated. Thus we believe that there could be other pathways than the method 
adopted by us to form glycine. These may be some of the reasons of the discrepancy seen 
in this Table. 

5. Concluding remarks 

In this paper, we presented the preliminary results on the chemical evolution inside a 
collapsing interstellar cloud. Our models are distinctly different from the other models 
used by several authors (Aikawa et al., 2005; Ceccarelli et al., 1996; Lim et al., 1999; 
Shalabiea and Greenberg, 1995; Shematovich et al., 1997). We do not include the heating 
and cooling processes while determining the dynamics of the cloud, and thus our model 
is not totally self-consistent. While in the outer edge of a diffused cloud the heating and 
cooling time scales are comparable to the infall time scales and should have been included, 
deep inside, the infall time scale is much shorter and cooling can be ignored. On the other 
hand, the outermost shell is also of very low density and thus the heating is low. The 
cooling is also negligible as the reactions rates are low. Thus we believe that even if the 
heating and cooling were included, the result would not have differed significantly. 

Unlike the previous study Chakrabarti and Chakrabarti (2000a, 2000b) we have incor- 
porated the grain chemistry of H2 formation self-consistently. This major improvement 
gave the most realistic abundances of H 2 molecules in the grain and the gas phases. 
We find that our computed average abundances generally agree with the observed abun- 
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dances (see, e.g., Allen and Knapp, 1978; Friberg et al., 1988; Irvine and Hjalmarson, 
1983; Matthews et al., 1985; Ohishi et al, 1995; Tolle et al., 1981) except when the 
molecules are complexes. We always seem to underestimate them as compared to the 
observed values or 'upper limits'. It is not clear at this moment how close the results of 
the bio-molecules are in comparison to the actual values, since the reaction rates we used 
are similar to the neutral-neutral rates which need not be accurate. If, for instance, the 
reaction rate at each step of glycine formation were higher by a factor of ten, the resulting 
glycine abundance would have been 10 3 times higher, which is closer to observed claims. 
Given that the pathways to produce glycine may itself be different, perhaps one needs 
to look into laboratory experiments which simulate interstellar clouds for guidance (see, 
e.g., Elsila, 2007). 

In future, we plan to improve the hydrodynamic model by including the angular motion 
and shock formation in the flow. We expect that jets and outflows would form and a part 
of this will fall back on the disk and the matter would be recycled. We also plan to 
improve the grain chemistry to include the formation of CH3OH, CO, NH3, OH on 
the grains themselves. Thus we anticipate that the chemical abundance will be strongly 
affected by such recycling of matter and such incorporation of newer species on the grain 
surfaces. 

A. Das acknowledge the support from an ISRO project. 
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SOLUTION PROCEDURE 

The equations are solved on a spherical grid extending from r in to r out composed of N 
equal logarithmically spaced grids along the radial direction. The code is customized to 
take care of the collapsing spherical hydrodynamical flow which is strictly one dimensional. 
That is, no back flow is possible. Thus, to solve Eqs. (1-2) we use the first order upwind 
differencing method. After appropriately splitting, the Eqs. (1-2) become, 



Pi +1 = Pi- r 2 , dt r , {pi+iVri+ir- +1 - Plv r irj), (Al) 

r i \ r i+l ~ r i) 



ri{r i+ i - r i} 

A+l V r{+1 r i+l/ Pi+1 ~ Pi v n r i/Pi) 



{^i+l ~ r i) 

Here, % denotes the index for the radial grid and j denotes the index for the time. To avoid 
the instability in the code we chose the time step by using the Courant-Friedrichs-Lewy 
stability criterion, which gives, 

|u|A*/Ar < 1, (A3) 

i.e., 

At ~ Ar/|u|, 

where, \v \ is the magnitude of velocity, At is the time step, and Ar is the grid spacing along 
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the radial direction. We always advance the time step after ensuring that the Courant 
condition is satisfied. To be on the safer side, we chose time step dt = At/2. Even though 
we use self-gravitating flow, we do not solve Poisson equation to get the potential </>j(r) 
here, since we are dealing with a spherical flow. Potential at any point is computed as 
a sum of two terms, one coming from the cloud itself [(fidoud = —GMd OU( i(r)/r, where, 
M c i ou d(r) is the mass of the cloud Ti n < r < r out within the grid] and the other is due to 
the cloud core 4> core = —GM core /r, where, M core is the mass of the core within r < r in , 
the inner edge of the computational grid. The value of <f>(r) that is required during 
the simulation is dynamically calculated from the mass within r, i.e., = —G(M doud + 
M core )/r. Here, M doud at jth time step is calculated by adding contributions from each 
spherical shell, 

M do J = EiAmrfdripi, (A A) 

where, i — 1,2, N and the densities at the jth time step are used. We put r = r out to 

get the potential at the outer boundary. 
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